{
 "metadata": {
  "name": "",
  "signature": "sha256:d0efcf30734ad00340cac1ab3445d22defe778f63c88b32bc3f2e578a901eee8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Computing a sparse solution of a set of linear inequalities\n",
      "\n",
      "A derivative work by Judson Wilson, 5/11/2014.<br>\n",
      "Adapted from the CVX example of the same name, by Almir Mutapcic, 2/28/2006.\n",
      "\n",
      "Topic References:\n",
      "\n",
      "* Section 6.2, Boyd & Vandenberghe \"Convex Optimization\" <br>\n",
      "* \"Just relax: Convex programming methods for subset selection and sparse approximation\" by J. A. Tropp\n",
      "\n",
      "\n",
      "## Introduction\n",
      "\n",
      "\n",
      "We consider a set of linear inequalities \n",
      "$Ax \\preceq b$ \n",
      "which are feasible. We apply two heuristics to find a sparse point $x$ that satisfies these inequalities.\n",
      "\n",
      "The (standard) $\\ell_1$-norm heuristic for finding a sparse solution is:\n",
      "    \\begin{array}{ll}\n",
      "    \\mbox{minimize}   &  \\|x\\|_1 \\\\\n",
      "    \\mbox{subject to} & Ax \\preceq b.\n",
      "    \\end{array}\n",
      "\n",
      "The log-based heuristic is an iterative method for finding\n",
      "a sparse solution, by finding a local optimal point for the problem:\n",
      "    \\begin{array}{ll}\n",
      "    \\mbox{minimize}   &  \\sum_i \\log \\left( \\delta + \\left|x_i\\right| \\right) \\\\\n",
      "    \\mbox{subject to} & Ax \\preceq b,\n",
      "    \\end{array}\n",
      "where $\\delta$ is a small threshold value (which determines if a value is close to zero).\n",
      "We cannot solve this problem since it is a minimization of a concave\n",
      "function and thus it is not a convex problem. However, we can apply\n",
      "a heuristic in which we linearize the objective, solve, and re-iterate.\n",
      "This becomes a weighted $\\ell_1$-norm heuristic:\n",
      "    \\begin{array}{ll}\n",
      "    \\mbox{minimize}   &  \\sum_i W_i \\left|x_i\\right| \\\\\n",
      "    \\mbox{subject to} & Ax \\preceq b,\n",
      "    \\end{array}\n",
      "which in each iteration re-adjusts the weights $W_i$ based on the rule:\n",
      "    $$W_i = 1/(\\delta + \\left|x_i\\right|),$$\n",
      "where $\\delta$ is a small threshold value.\n",
      "\n",
      "This algorithm is described in papers:\n",
      "\n",
      "* \"An affine scaling methodology for best basis selection\"<br>\n",
      "  by B. D. Rao and K. Kreutz-Delgado\n",
      "* \"Portfolio optimization with linear and fixed transaction costs\"<br>\n",
      "  by M. S. Lobo, M. Fazel, and S. Boyd\n",
      "\n",
      "## Generate problem data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cvxpy as cvx\n",
      "import numpy as np\n",
      "\n",
      "# Fix random number generator so we can repeat the experiment.\n",
      "np.random.seed(1)\n",
      "\n",
      "# The threshold value below which we consider an element to be zero.\n",
      "delta = 1e-8\n",
      "\n",
      "# Problem dimensions (m inequalities in n-dimensional space).\n",
      "m = 100\n",
      "n = 50\n",
      "\n",
      "# Construct a feasible set of inequalities.\n",
      "# (This system is feasible for the x0 point.)\n",
      "A  = np.matrix(np.random.randn(m, n))\n",
      "x0 = np.matrix(np.random.randn(n, 1))\n",
      "b  = A*x0 + np.random.random((m, 1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## $\\ell_1$-norm heuristic"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create variable.\n",
      "x_l1 = cvx.Variable(shape=(n,1))\n",
      "\n",
      "# Create constraint.\n",
      "constraints = [A*x_l1 <= b]\n",
      "\n",
      "# Form objective.\n",
      "obj = cvx.Minimize(cvx.norm(x_l1, 1))\n",
      "\n",
      "# Form and solve problem.\n",
      "prob = cvx.Problem(obj, constraints)\n",
      "prob.solve()\n",
      "print \"status:\", prob.status\n",
      "\n",
      "# Number of nonzero elements in the solution (its cardinality or diversity).\n",
      "nnz_l1 = (np.absolute(x_l1.value) > delta).sum()\n",
      "print 'Found a feasible x in R^{} that has {} nonzeros.'.format(n, nnz_l1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "status: optimal\n",
        "Found a feasible x in R^50 that has 40 nonzeros.\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Iterative log heuristic"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Do 15 iterations, allocate variable to hold number of non-zeros\n",
      "# (cardinality of x) for each run.\n",
      "NUM_RUNS = 15\n",
      "nnzs_log = np.array(())\n",
      "\n",
      "# Store W as a positive parameter for simple modification of the problem.\n",
      "W = cvx.Parameter(shape=(n, 1), nonneg=True); \n",
      "x_log = cvx.Variable(shape=(n,1))\n",
      "\n",
      "# Initial weights.\n",
      "W.value = np.ones((n, 1));\n",
      "\n",
      "# Setup the problem.\n",
      "obj = cvx.Minimize( W.T*cvx.abs(x_log) ) # sum of elementwise product\n",
      "constraints = [A*x_log <= b]\n",
      "prob = cvx.Problem(obj, constraints)\n",
      "\n",
      "# Do the iterations of the problem, solving and updating W.\n",
      "for k in range(1, NUM_RUNS+1):\n",
      "    # Solve problem.\n",
      "    # The ECOS solver has known numerical issues with this problem\n",
      "    # so force a different solver.\n",
      "    prob.solve(solver=cvx.CVXOPT)\n",
      "    \n",
      "    # Check for error.\n",
      "    if prob.status != cvx.OPTIMAL:\n",
      "        raise Exception(\"Solver did not converge!\")\n",
      "\n",
      "    # Display new number of nonzeros in the solution vector.\n",
      "    nnz = (np.absolute(x_log.value) > delta).sum()\n",
      "    nnzs_log = np.append(nnzs_log, nnz);\n",
      "    print ('Iteration {}: Found a feasible x in R^{}' + \\\n",
      "           ' with {} nonzeros...').format(k, n, nnz)\n",
      "\n",
      "    # Adjust the weights elementwise and re-iterate\n",
      "    W.value = np.ones((n, 1))  \\\n",
      "              /(delta*np.ones((n, 1)) + np.absolute(x_log.value))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Iteration 1: Found a feasible x in R^50 with 48 nonzeros...\n",
        "Iteration 2: Found a feasible x in R^50 with 36 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 3: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 4: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 5: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 6: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 7: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 8: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 9: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 10: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 11: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 12: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 13: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 14: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Iteration 15: Found a feasible x in R^50 with 33 nonzeros..."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Result plots\n",
      "\n",
      "The following code plots the result of the $\\ell_1$-norm heuristic, as well as the result for each iteration of the log heuristic."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "\n",
      "# Show plot inline in ipython.\n",
      "%matplotlib inline\n",
      "\n",
      "# Plot properties.\n",
      "plt.rc('text', usetex=True)\n",
      "plt.rc('font', family='serif')\n",
      "plt.figure(figsize=(6,6))\n",
      "\n",
      "# Plot the two data series.\n",
      "plt.plot(range(1,1+NUM_RUNS), nnzs_log, label='log heuristic')\n",
      "plt.plot((1, NUM_RUNS), (nnz_l1, nnz_l1), linestyle='--', label='l1-norm heuristic')\n",
      "\n",
      "# Format and show plot.\n",
      "plt.xlabel('iteration', fontsize=16)\n",
      "plt.ylabel('number of non-zeros (cardinality)', fontsize=16)\n",
      "plt.ylim(0,n)\n",
      "plt.xlim(1,NUM_RUNS)\n",
      "plt.legend(loc='lower right')\n",
      "\n",
      "plt.tight_layout()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9sHPlZ+PG3feVn28R2VK7XVr1kk4NSBJfYvh4IKsI5\nNqIFtWoSHz8ESHBO8kXiH9TEuUqoKRK9/EBISKA7O8evfxCJnQooAi6xjxREUTknvkNVUXUXpwel\nvZbGdg5KD0Fvvn98drzr9azt2Z31zo7fL8na3dmZ2Wdbx899PvPM8wFJkiRJkiRJkiRJkiRJkiRJ\nktrqXPlxrGrbYWCoZpskaRvrbsNnjgEvArfKr/vLj7PlxwNbHpEkKXfalaAeAJ4tvx4FlsrPF4BD\nbYhJkpQz7UhQfYTpvJPl1z3AYtX7u7Y8IklS7ryhDZ95sfw4TEhUAF3rHfDggw9GL7zwQkuDkiS1\nzQvA/tqNWz2CGiMURADcAUrAMmFUBdBb3r7KCy+8QBRFW/Lz0Y9+dMs+ayt//F6d81PE7+T36ryf\nrfxewINJCWOrE9QCMFN+vgt4DrhESFQAe4BrWxyTJCmHtnqKb5bKCOprwPPl54OE6b7lqm2SpG2s\nHdegriRsi69LzSa8t6UOHjzY7hBawu/VOYr4ncDv1Wny8L3WLU7Ikag8TylJKpiuri5IyEftKDOX\nJGlDJihJUi6ZoCRJuWSCkiTlkglKkpRLJihJUi6ZoCRJuWSCkiTlkglKkpRLJihJUi51TIL6+tfb\nHYEkaSt1TIK65iIckrStdEyC+uQn2x2BJGkrdUw383vvjfjSl6C7Y1KqJGkzOr6beV8fzM21OwpJ\n0lbpmAT1Uz/lNJ8kbScmKElSLnXMNaj/+7+It74VbtyAd76z3eFIkrLS8deg7rkHfuIn4C//st2R\nSJK2QsckKHCaT5K2k46Z4ouiiFdfhXe8A770JXjTm9odkiQpCx0/xQewYwc8/LBdJSRpO+ioBAVO\n80nSdtFRU3wAt2/DD/4gfPnLdpWQpCIoxBQfwJ498Ja3wD/9U7sjkSS1UsclKHCaT5K2AxOUJCmX\nOjJBPfwwvPIKvPxyuyORJLVKRyaoe+6B973PUZQkFVlHJihwmk+Siq7jysxj//mf8Pa3w7//O7z5\nzW2KSpLUtMKUmcfe/Gb4oR+Cq1fbHYkkqRU6NkGB03ySVGQdO8UHoYrvoYdCV4l77mlDVJKkphVu\nig/g/vvhrW+Fz3ym3ZFIkrLW0QkKnOaTpKIyQUmScqnjE9R73gP/8R+hy7kkqTg6PkF1d8P73+8o\nSpKKpuMTFDjNJ0lF1NFl5rH/+i9429vgi18My8JLkjpHIcvMY296E/zwD8Mzz7Q7EklSVgqRoMBp\nPkkqmkJM8QH867/CwEBYJ8quEpLUOQo9xQfwzneG7ub/+I/tjkSSlIXCJChwmk+SisQEJUnKpUIl\nqMFBWFqCW7faHYkkqVmFSlB2lZCk4ihUggKn+SSpKApTZh77+tfhvvvg3/4Ndu5scVSSpKYVvsw8\n9sY3wnvfC3/zN+2ORJLUjMIlKHCaT5KKoHBTfBCaxj74IHzlK/CGN7QwKklS07bNFB/AO94B998P\nn/50uyORJDWqkAkKnOaTpE5ngpIk5VJhE1R/P7z6Krz4YrsjkSQ1Im2RxGFgECiVXy8AzwGfyDKo\nBKmKJGLHjsG73gW/9mstiEiSlIl6RRKbTVAngceBnjrvLwHjwNONBLcJDSWoT34Sfvu34W//tgUR\nSZIy0WiC2gnMEkZKl8rPl2v2KQH9hAT2EvBok7EmaShB/fd/w1vfCi+/DL29LYhKktS0RhPUZeAJ\nYH6Tn3OEMAV4Ok1wm9BQgoJQLPGzPws/8zMZRyRJykSjCaqHtSOmjewE7qY8ZiMNJ6jJSbh+Hf7k\nT7INSJKUjUZv1K1NTo9s4rM2m5xOVj0/DAwBY5s8dtN+8idDX77//d+szyxJaqW0ZebTwP4MPvcQ\nMFx+3l9+nC0/Hsjg/Cve9jYoleAf/iHLs0qSWi1tguoBzgNXgcea+Nzq+bpHCVWAEIoxDjVx3kTe\ntCtJnSdtgjoNjACj5WPngCdJN6o6QGW0BOGa1WLV610pY9qQCUqSOk/aBHW+/LgMTAJnCSOgm4RR\n1WauUfUlbGtpV/UDB0LJ+ec/38pPkSRlKW2CugTsIBQ4LBLK0BeB44RR1SAhUe2oc3zt6AlCsouT\nVi9wJ2VMG+rqCsUSjqIkqXOkHbm8XvV8muR7pPoJXSWSbtg9XH7cBRyjUrU3CFwkJL5rwPM1x0X8\naNWr3cAe+OiPfpQzB8+s+ZAz18/wsU99bM32d97+KC//0eb3T3t+93d/93d/99/E/n/0MfhC1cZP\nAU20Ooq9Tpjme4LkcvKdwBQh4SRN5cXGgFPAUUIyGiMUSJQIiapWw/dBxb7xjdBV4vZt6FsvMknS\nlmq2F1/sLOt3iegBbgBPARdSnns9TScogA98AEZH4ed+LoOIJEmZyGpF3Ws1r/sJlXzxtallYC/Z\nJqfMWM0nSZ0j7QjqKqHMvNYhwrTfQ01HlCyTEdSXvwzvfjd89avwLd+SQVSSpKY1M4LaSRgd7Sy/\n3pHwsxMYyCLQVrrvPnjgAfj7v293JJKkjbxhg/eHCBV51d0d6jWPnc4kohaLp/ke2cwdW5Kktkkz\nxTdBSFjnao5bJCStmQzjqpXJFB/A88/DkSNhKfiult4eLEnajKyq+M4RRlRbLbMEFUVw//3wzDPw\nvd+bySklSU3IqopvveT0oZTnagu7SkhSZ9hoBBU3gY07O9S7ctNFuPfpgSyCSpDZCArgr/8aPv5x\niyUkKQ8aneJbAr5GJfG8vs6+EXBPI8FtQqYJ6rXX4N57YWEBdmXeO12SlEa9BLWZKr7qqr1lQnui\npZr9+ggjqI7w7d8eqvj+6q/g53++3dFIkpJslKBu1ryeYm038lg7iicaFpebm6AkKZ+yLLQ+DFzJ\n8HzVMp3iA/jKV+Bd7wqP3/qtmZ5akpRCo1N8B1i9PHs9uwiNZFuVoDJ3773wPd8Df/d3cCjzReYl\nSc3aKEHdSHGubIc4WyCe5jNBSVL+bDTFt0goitjMVOBTwL6mI0qW+RQfwD//M3zwg3Drll0lJKld\nGp3iW68oolZHFUkAfP/3wze/CZ/7HHzf97U7GklStY06SRxPca56TWRzq6vLNaIkKa/StjqK7SZ0\nmYh/DgOXM4ppS5mgJCmf0l55OUyY9ksyQ/JihlloyTUogP/5n1DR9+KL8Ja3tOQjJEnryKpZ7OPA\neWCUcBPvUeAEsAAcaS7E9vi2b4OhodBVQpKUH2kTVA9wmrA44RyhDH2SMHK6mG1oW8dpPknKn7QJ\naqHq+QyVlXYX6IAl3+t53/tgZiZM90mS8mGjMvNad4GrhJtyR6ncyDtIaBjbkb7ru+Dd74ZPfQpG\nWnUVTZKUStoR1BjhQtbd8s9pwhTfsfJjx3KaT5LyJYv+CSXCtanazudZalkVX+yznw0r7d6+bVcJ\nSdpKWVXxJVkgJKd6q+12hLiTxGc/2944JElBMzfqVv8cAM5lEVC7dHWFtaFOnICvfrXd0UiS0hZJ\n7CEURvQkvNdx3cxrfexjIVE9/DD8+Z/DD/xAuyOSpO0r7dWWufLjBGuXfT9Lh3Uzr+dP/xR+9Vfh\n6afhAx/Yso+VpG2p0W7mtfqBXkIFX63e9GHl00//NJRK8KEPwb/8C4yPWzghSVst7TWoWerfkLtQ\nZ3tHes974DOfgelp+MVfhNdea3dEkrS9pE1QZwltjj5MqNrbTaVI4myWgeXB298eloR/7TX4sR+D\nV15pd0SStH2knbh6fZ33IuCeJmJZz5Zeg1r74fAbvwF/8AeheGL//raFIkmFU+8aVNoEtQQ8Vue4\nwhRJ1DM1Bb/yKzAxEa5PSZKal1WRxBPAlTrvtT+DtNjRo6F44oMfDMUTH/mIxROS1CpZ/nk9TP3k\n1axcjKBiX/pSSFL79sHv/z58x3e0OyJJ6lyNTvHFV1ueLz/Wa2fUBTwFPNBIcJuQqwQF8I1vwC/9\nEiwswJ/9Gdx3X7sjkqTO1GiCWgK+RiXxbMsiiXqiCH7zN2FyMiSp/v52RyRJnafRBNUPLFO5x2mR\nsMx7bReJPsIIqtBFEvVcuRJ6+D35JBzpyIXvJal9Gi2SqF1CY4pws26S8fRhFcPhw5Xiic99Dn79\n1y2ekKRmZfln9BLwaIbnq5brEVTslVdCkrr/fvjDP4Tv/M52RyRJ+dfoFN9ZNlc+3ktYbXdbXYNK\n8tpr8Nhj8PnPh+tSb397uyOSpHxrNEGtVxSRJIsFEJN0TIKCUDxx7hz83u/BJz4BDz3U7ogkKb8a\nXVF3obxP/DNCpWFsX3nbPuAioXhChOtPp0/D7/4uvP/9cOlSuyOSpM6z0QjqJHCh6vWL1L/X6Rng\nx7MIKkFHjaCqvfBCWFPqF34BzpyB7laNMSWpQ2XVi2+RUHr+hZrtOwkr7W7LMvONfOUroXfffffB\nH/8xvPGN7Y5IkvKj0Sm+WhcJ035PEpbceIxQSHEDmGkuxOK691549tmQmN77XvjiF9sdkSTlXyNl\n5hOEir1q08Bo8+HU1dEjqFgUwYUL8Du/E4onHn643RFJUvtlNcUX6wFKhEKJOUK3iVYqRIKK/cVf\nwC//Mnz3d7c7Eklqv09/OpsEdRI4RqjYe36DfbNUqAQFYZrv5ZfbHYUktd+P/Eg2CeoW4cbdo8B8\n82FtWuESlCQpyKpIYgIYJjk51V6XkiSpYWlHUHuA4+XHa1SuP3URGskOZhpdhSMoSSqorIokXA9K\nkpSpRpfbSHI06USE+6EkScpE2gR1AbhS573eJmORJGlFlutBDVF/McNmOcUnSQWV5RQfwG7Czbqx\nvcAksKvB80mStEraBHWYUK2XxF58kqTMpJ3imyMkojngNPAEYdR0itDl/NVMo6twik+SCiqrKb4e\nQmICOEToYv4FQtK6CDzacISSJFVJ20lioer5DCFJxdsHMolIkiTSj6DuAlcJN+WOEkZQEDpI9GUY\nlyRpm0uboMYIRRLLhGR1Grhcfu98hnFJkra5LO6DKhGWfG9ld3OLJCSpoLLqZg5hqfcPVb0+QGge\nK0lSZtKOoJ4iLFi4xOqbcqeAzwC/lVFctRxBSVJBZdXNfI7QLLYXuFm1vYdQybeZQokjhAR3FDhR\n3naYcF2rRChXr2WCkqSCynKK7zarkxNsvsR8iErPvhJherC//F7cx+9AAzFJkgombRXfLPAicI7K\nPVEDwONsrtXRLJVE1EcorDgHPFPetkC4t2orl5OXJOVQ2gQ1Thj5TNZsv8nml3zfSbiO9UTV68Wq\n9204K0lqeMHCEmFqrg94jnQjnruEdaWuUpkq3PBa2JkzZ1aeHzx4kIMHD6b4SElSXly/fp3r169v\nuN9GiaGRNZ7WO6af0IVinrAC7x3CiOla+ZgjhJL1CzXHWSQhSQXVaJHECJWpuI3sJHSV6FlnnyEq\nlX49wC3gEmFEBiE5Xdvk50mSCmyjKb5xQhHDIiGR3CAUMsTXjPoIyWWEUNwwTv0l4SFcuxotH7ME\nfKK8fZCQvJaB59N+CUlS8Wz2Pqg4+QzVeX+y/P7dLIJK4BSfJBVUVjfq9hBGOyXCKOomq5fgaBUT\nlCQVVFYJql1MUJJUUFl2kpAkqeVMUJKkXDJBSZJyKW2COgDsB3aXX48R7n36cIYxSZKUOkFdJFTu\nHQVOAhOEar5u4MlsQ5MkbWdpe/EtEJLTbcLaUMtU1nSayzAuSdI2l3YE1UNIThD66l2uem9x7e6S\nJDUm7QiqD9gBDJdfx33zDrN2EUNJkhqWNkGdJiSiEmGBwhlCkhrCJq+SpAxl0UkiXrI9XkajFewk\nIUkFVa+TRCMLFkKY5jtUfj4DvNrgeSRJStTIjbpPEar3pss/y1hiLknKWNopvrOEkdMlKtN5JeA4\nYQn3x7MLbRWn+CSpoLLqZj5HWG4j7XvNMkFJUkFl1c18vXudvA9KkpSZtEUS88AzhBZH8UKFewlT\nfN4HJUnKTCNl5lOEG3OrTQOjzYdTl1N8klRQWa+oWyLc/9RHuEH39vq7N80EJUkFlVWCOgkcIzSM\nfb75sDbNBCVJBZVVkcSJ8kmy6EAhSVJdaRPUBKFRbFJLo7Hmw5EkKUg7EtpDqNjbQ7j2FK8J1UUo\nnvA+KElSKlldg3p9nfci4J6U59ssE5QkFVSWzWKPJp2I0AZJkqRMpE1QF4Ardd7razIWSZJWNFqN\nt5twL9SzhOtR3gclSWpIVmXmELqWLxAq+gAGCMUSOxoNTpKkWmkTVLwW1Chwt7xtmnDz7lSGcUmS\ntrm016BKwEj5+bGq7TeBhzKJSJIk0o+g6hVCHGg2EEmSqqVNULPAc8BjhGT1CKE/3w3gcrahSZK2\ns0aq+CZY29ZohsrUXytYxSdJBdXK5TbmaP1ihSYoSSqorBPUVjNBSVJBZXUf1HPA/pptV4GXgBcb\nikySpASNNIuNgEng/1VtHyAkr0Zu/N0MR1CSVFBZjaBuAr9FWHLjDqGKD0IVX6uvQ0mStpG0CWoR\nGCes+7RMqN57EthZfk+SpEykneJbJDSKfbX8+hzhPqhlwtTfrswiW80pPkkqqKym+HqALxBGTRBG\nU/sI3cx7Gw9PkqTV0vbiS0poC4QiiSPNhyNJUtBM1d0lwnRfbLq5UCRJqmg0QfUQln4vZRiLJEkr\nWnXfkiRJTTFBSZJyqdEEFa+qeyPDWCRJWpFls9ghwnpRreB9UJJUUFl3M99d87qX0J+vVcu+m6Ak\nqaDqJai090HtIUzr9SS8ZwaRJGUm7Qhqrvw4ASzVvHeW0FWiFRxBSVJBZTWC6idM591NeM9WR5Kk\nzKSt4psltDVKstBkLJIkrUg7xTcETAEfJ6z/FCcliyQkSQ3Jqorv9XXei4B7Up5vs0xQklRQWV2D\nugs8lnQiQpGEJEmZSJugngCutCIQSZKqNXqj7g7gUPn5DJUVdlvFKT5JKqisVtQFeIrQi2+6/LNM\nZYVdSZIykXYEdZYwcroEzJe3lYDjwFXg8exCW8URlCQVVFZVfHPAYAPvNcsEJUkFldUU32KD70mS\nlEraKr554BlCL774Jt29hCm+mxnGJUna5hqp4psCDtdsmyYsYNgqTvFJUkFlvR5UidA4tg+4BtxO\ncexY+XEvcLr8/DChGrAEXEw4xgQlSQWVVSeJ2AJrm8NuZkXdIcJ9U7eBy+XX8bWrWUKCOkClQlCS\ntE01ch8UhBV1q38OsLlWRyUqN/gulF8/Shg9xdsOJRwnSdpmtnpF3erpu37C/VQDwJ2q7btSxiRJ\nKqC0CWqKMMqpt6LuZvUTEl08ldfotTBJUkG1a0XdISpdJ5YJxRbxOe4kHXDmzJmV5wcPHuTgwYMp\nPk6SlBfXr1/n+vXrG+6XduRyjdDR/NmE9zZTJAFwjLC4YXzMIqEDxUXgZPkznq85xio+SSqorDpJ\nnCXc8/Rh4BHSF0kcKu/3EiExRVSm+YYIo6na5CRJ2oZcUVeS1FauqCtJ6iiuqCtJyqVOKe92ik+S\nCirLFXUlSWo5E5QkKZc2SlBjwHOEUnJJkrbMRkUSx8uPnXKtSpJUEBuNoOIuD/F6T5fW2fexTCKS\nJImNR1BdhNLya4TmsHuB/XX2Ow08nWl0kqRta6Opux5C1/E9mziXnSQkSak1u+R7vLz7WWA84bie\n8nv7Gg9xXSYoSSqoZlsd3Sw/zrC5juWSJDWl0eq8nYTu4xCS1qvZhFOXIyhJKqgsO0k8RSiYmC7/\nLANPNhOcJEm10o6gzhLWdLpEZR2nEuF+qatUVsnNmiMoSSqoZoskYnOE+6LSvtcsE5QkFVRWU3yL\nDb4nSVIqadeDmgeeASaAhfK2vYQpvpv1DpIkKa1GqvimgMM126aB0ebDqcspPkkqqKyuQcVKVG7e\nvUalV1+rmKAkqaCyTlBbzQQlSQXlirqSpI5igpIk5ZIJSpKUS2kT1EngRZLXhJIkKTNpE9QJwoWs\nTimukCR1qLQJagIYptKHr9pY8+FIkhSkHQntIXSN2EO4/2mO0M28i3ADr734JEmpZHUf1OvrvOeS\n75Kk1JpdUbfa0aQTEZbikCQpE2kT1AXgSp33+pqMRZKkFY1W4+0m9ON7lnA9yl58kqSGZNnq6Cph\nqY2J8usBQrHEjkaDkySpVtoE9RSham8UuFveNg0cI1TxSZKUibTXoErASPn5sartN4GHMolIkiTS\nj6DqFUIcaDYQSZKqpU1Qs8BzwGOEZPUIoT/fDeBytqFJkrazRqr4Jljb1miGytRfK1jFJ0kF1col\n3+cI16BayQQlSQXViiXfDwBLwBeaOMdmmaAkqaCyvA/qKUJPvhuE+6G+CXy4meAkSaqVtsz8LKG8\nfJKQnAD2Ah8pP/+tjOKSJG1zaaf4FgmdI2pbG/UQCiVcbkOSlEpWU3wLJPfdWyYkL0mSMtHIirpJ\n15vGqEz5SZLUtI2m+F4iLERYvX+pvG2h6vUSoav5qy2IEZzik6TCarTM/HVCM9h6+/UQ7oXqKe/7\nQOMhrssEJUkF1WiCugYMtyKglExQklRQjRZJpElOT6QJSJKk9TTSSWIIOJRwnjFgV9MRJXMEJUkF\nVW8ElfZG3ZPAuTrvmUEkSZlJW2b+ODBePq7259lsQ5MkbWdpp/heInQxTyonPwDMNx1RMqf4JKmg\nsuokcR64COxIeO9s+rAkSUqWdgS1k9DqaCer2xvtKm+7J7vQVnEEJUkFlVWRxBThptwrCe8dTh+W\nJEnJ0iaoQ4Ru5knXmp5qPhxJkoK016DmCX33kkw1GYskSSvSjqA+Tmh/dA6Yq9reRSiSeCijuCRJ\n21zaIonX13kvwiIJSVJKWRVJ3AUeKz9W68Eyc0lShtImqEmSK/jAVkeSpAylLZIYb0kUkiTVaKSb\neZI9hKIJu5lLklLJ6hpUUpFElHRiSZKakTZBQZjmixPSLkLz2AHSLVh4jtXThYcJrZNKhF5/kqRt\nLm2CmgcuJGzvB0Y3eY5jhIQUJ6j+8uMsIUG1siu6JKlDpC2SGKiz/SaVRLORSWCh6vUole4UC6xd\nrVeStA2lHUElLbMBMEwY/TSih0pXdGhdoYUkqYOkTVDL67zXTAm6RRaSpFUaKZI4zuqEskiYmrvZ\nYAzLQF/5eS9wJ2mnM2fOrDw/ePAgBw8ebPDjJEntdP36da5fv77hfmlHLrXVd426CoyUnx8ABgnV\neycJzWifr9nf+6AkqaCyWvJ9veQ0tMlzHCEkpMfKr+OKvSHCaKo2OUmStqFGr/3srnndS6jOa9Vy\nG46gJKmgsuoksQe4Qai8q2UGkSRlJu0IKl6kcIK1K+ueBfY1HVEyR1CSVFBZjaD6CdN5tetBUd4u\nSVIm0hZJzFK/m8RCne2SJKWWdopvCJgCPk647ylOShZJSJIaUm+KL22CSlpuIxYB96Q832aZoCSp\noLK6BnWXcP9SUmI7mz4sSZKSpU1QTwBXWhGIJEnVOqVJq1N8klRQWbU6kiRpS5igJEm5ZIKSJOWS\nCUqSlEsmKElSLpmgJEm5ZIKSJOWSCUqSlEsmKElSLpmgJEm5ZIKSJOWSCUqSlEsmKElSLpmgJEm5\nZIKSJOWSCUqSlEsmKElSLpmgJEm59IZ2ByCps/X19bG0tNTuMNQBent7WVxc3PT+a9aAz6koiqJ2\nxyApQVdXF/771GbU+13p6uqChHzkFJ8kKZdMUJKkXDJBSZJyyQQlqZDOnz9PX18fFy9ebPpc4+Pj\ndHd3Z3Ku9YyOjm7qM27fvs3t27cbOraTWMUnqZBOnTrFzMxMfAG+KefOnWN+fj6Tc63nIx/5CD09\nPRvuNzU1xd69e9mzZ0/qYzuJIyhJyon9+/eze/fudfdZXl7m8uXLDR3baUxQkraNyclJTpw4wejo\nKFeuXEncfvz4cUZGRpidnV1z/NLSEqOjo2umDqenp1eOv3DhAhCmGLu7u7ly5Qq3b99meHiY7u7u\nVe9dvHiRixcvsm/fPi5cuEBvb+/K8QDHjx/n9OnTnDhxghMnTgBw+fJlFhYWuHTpEqdPn145X19f\n36pjJycnGR0d5cSJE4Wb+subSFI+bfTvE7L5acTw8HB08eLFKIqi6MaNG1Fvb+/Ke11dXdHt27ej\npaWlqKura9X22dnZaHl5edW5Dh06FA0ODkZRFEU3b95cOebWrVvR3r17V/br7e2N5ufnoyiKooGB\ngejKlStRFEXR8vLyqs8ZGBiIRkZGouXl5WhkZCSKoig6evRodOHChSiKoujatWvR+Ph4FEVRtLCw\nEJ04cWLl2L17966cN1Z97I0bN1Y+a3JyMtq3b1+K/9Vap97vCpB4I53XoCS1VF7u4b106RJ79+5d\ned3f38/U1BT9/f2USqWV7aVSib6+Pnbu3Lnq+K6uLh599FEADhw4AMDdu3eZnp4GWBnNPPTQQyws\nLLB///5Vx0cJ/0OUSiV27tzJM888s2a/UqnE+fPnmZmZ4dChQ5w6dSrVdx0YGABgbGyM4eHhTR+b\nJyYoSdtCV1cXvb29K6+jKKKrq4vBwUEWFxe5e/cuURSxtLS0JrnEaosQurq6uHPnDj09PZw9e3bd\nz19YWFizLU4iteeMH2/dusXMzAxTU1MMDw/z0ksvrdl/fn5+JWHW+66dem3Ka1CSCiuKopURyaOP\nPsrc3NzKe/Pz8xw5coQoihgaGuLIkSMcO3aMmZmZDc9VvW14eJibN2+ulH1PT0+vXL/q6+vjzp07\nAKs+O1adRGo/49q1a0xMTDA2NsbVq1dXJcdSqcSdO3eYmZlZSXzVxx46dIiZmRnu3r0LVEZ3ao2W\nz41Kakxe/31OTU1Fvb290eDg4Mo1ocnJyej48ePR0aNHV67hxNdruru7o66urqirqyuanp5OPNfI\nyEi0sLAQnTp1Kuru7l65LnT+/PloYGAgGh4eXrkOFEVRNDMzEw0MDETj4+PR5ORk1NXVFY2MjETT\n09Mr57t58+aaeG/evLkS6/j4eHT8+PGVa2nx9+jt7Y1GR0frftfz589He/fujQYGBla2tVu93xXq\nXIOyWaydPgJBAAAHhElEQVSkpnR6s9h4xPPkk08CMDs7y7lz57h69WqbIysem8VKUgqXL1+mv79/\n5XVvb++qYgq1jyMoSU3p9BHU/Pw84+PjlEolenp6WF5e5vz58+zYsaPdoRVO2hGUCUpSUzo9QWnr\nOMUnSSoEE5QkKZdMUJKkXDJBSZJyyQQlScolE5SkQqq3ou7y8jLDw8PMz8+3KbKNuYJvYLNYSYWU\ntKLuhQsXmJiYYGFhoeWr4zbDFXwDR1CSto2TJ0+udAQv2vLojcj7Cr4mKEktdeb6Gbo+1rXm58z1\nM5vev96+rTI5OUl3dzcnTpxgZGRkzVRh0sq89VbJrT3PlStXOH78OIODg4yOjq4bhyv4doYt6LMr\nqRF5/vdZvaJutXg13Y2OTVpBt97KvFGUvEpu9Xnijuazs7NRFIXVd2dmZhI/v4gr+Nb7XaFON3NH\nUJK2vXhE0NfXx9NPP72yPWkF3aSVeeNVdaMoSlwlNz5PvEDhI488AoR1nWqLD2KbWcH39OnTKyv4\n1ooaXMF3cHCQiYmJplbwvXbt2qaPXY9FEpK2vVOnTq1MaVU3iU1aQTdpZd7q95NWya0+T/Xzvr6+\ndePa7iv4OoKSVFhRwiq4saWlpVWvd+zYsSo5JR0bRVHdlXnj99dbJTdNfPU+3xV88ydx3lJS++X1\n32e9VWaPHj0adXd3r6x0u96x9VbQTVqZt/qYeqvkHjlyJOru7o4uXry4sqru4OBgtLCwkOrzO3UF\n33q/K7iirqRWcLkNbZbLbUiSCsEEJUnKJROUJCmXTFCSpFwyQUmScskbdSU1pbe3N9edwZUftfdh\nbSQvv1WHgWWgBCR1GbTMXJIKKs9l5v3lx9ny44F6O26F69evt/PjW8bv1TmK+J3A79Vp8vC98pCg\nRoG458gCcKiNseTi/5RW8Ht1jiJ+J/B7dZo8fK88JKgeYLHq9a52BSJJyo88JCjIz7UwSVJO5CEx\nnAWuEa5BHQH2ABdq9nkeeHCL45IkbY0XgP3tDiLJAWCs/PwkOQ1SkrS17ml3AMArwADQC3wH8Fft\nDUeSJG2lk+0OQNvSuZrXh4EhKrMmnar2e42Vf9Zf5jb/ar9XzL8fbVaUX7Akh4Cr7Q4iY/2EP3ad\n/oeuVlH+gAMcA6rXDI//P4Pw/dp6z2MTar/XEOHaOcDl8utOVPu9Ym37+5GXKr52GwJmCF0sSnTu\nL1g9RWzDcRq4QrhNoVP/0NU6QLgXcLb82Onfa5LwPWK5uuexCbXfq0TluyyUX3ei2u8Va9vfDxNU\nUJRfsCQHqHTpKIojwHPl5xeA+TbGkrV4iqVEsb4XFPeex4tUWrT1U/ndLIK2/v0wQQVF/gXra3cA\nLTBI+ON2gGLNjc8Dtwl/xBc32LdT5eHWllbpB24Qbospirb+/TBBrVa0X7Aijp5iX6Mywji83o4d\npIdwDWCM8B9Me9bfveMsU/mD1wvcaWMsrTAEPN7uIDLU9r8fJqjVivYLViL88T5G+MPQ6dc0YncI\nIw0If/QeamMsWRoDJgjX1o4SpjKL5BKV6fM9hBv0i+IYlQYDRbmG3fa/HyaoiiL+gl0p/0TATopT\nLDFN5Q9dD/BPbYwla6+WH2cJybeTHSFMxz5Wfh2PeIcI361TZypqv9chQvXvS4Sp2U79d1b7vYr6\n96PjHCL8YsW/YI+0Nxxtwhjhv+6eaHcgGTtJMcvnJUmSJEmSJEmSJEmSJEmSJEmSJGWvn9B5u3aF\n53bdxJ30ufVilArBThJSsmVCv79a44S76rda0ufWi1GStA29DuzYRp8rtY0jKGnzJsqPW71kRLs+\nV5KUMxOEnoyvU+nLeI7Qq/F1wvLXV1l9XegcMFfe/lTN9tfLP3uAU8AtKutY9QBT5eNeKp+jevmQ\nuAlp7ec+VRVj7fWpc+X95srHx04RrlnFx1wrP5+jeEt7SFJhDbE6QUFo4Jo01TbB6kUuLxMSROxw\n+bizhETwEvBi+b1rNfueLe9brd7nJsU4BTxT9XqO1QkzjiVuiryzvE9RF0hUB3OKT0q2lLAtaYqt\nRGUdp9hlQof83TXH7SGsY3W0/BN/TvUyBnPlx91V2+pN7dXGGK/fUx3LJcJSMnFyi8/1ceBZ4G55\n/56az5Ta7g3tDkDqcP3lxxPAaPl5H2EarwR8oWrfeIG++apto4RRzBFCUhssb+9pIpaFqm3xZx0C\nPlG1/WbV86RkLLWdIyipcf1UljA/BYyUfwaBBwgjlGpzrHUEuEFISGfZ3PpW/XW299XZnmRh412k\n9jJBSZtXe51mksq1p4Ga946wtvCgdoXcEmE68CzwNGG0tWuTn5tkpvy4t+Yzqt+TOoYJSkoWj0aq\nr//ECWaYMGV2izCFNkOonDtQfr8fOE243lSt9lpSPI0XT+uVCCMxWJ1kkj43KcYFYBo4XnXsccI1\npldZrZe1krZJknLkCKHS7puEarsPVb13mVAF9xyrq+risu9FQhVd/N6xqnPNsXYp97PlY14iFDTs\nIFT13WF1dV7t564X41kqZebVU4bVxzxHqAIcS9gmSZIkSZIkSZIkSZIkSZIkSZIkSZIkSepU/x84\nYYxkXT7qswAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x106bb8c50>"
       ]
      }
     ],
     "prompt_number": 4
    }
   ],
   "metadata": {}
  }
 ]
}
